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ABSTRACT 


A classic problem in the study of the physics of sound 
in fluids is that of the finite amplitude plane wave propaga- 
mcg in a viscous, unbounded medium. Though the solution to 
fees problem 1s well known, for a given range of parameters, 
it would be desirable to develop techniques for solution of 
the problem which are not similarly limited. Through the ap- 
Peteation of the technique of parametric differentiation, 
this goal is realizable. This extension method transforms 
the governing non-linear differential equation to a linear 
equation in what may be termed parameter space; the equation 
is solved and a guadrature recovers the dependent variable 
sSOlution. Parametric differentiation is conceptually straight- 
forward and has been applied in the past to a wide variety of 
non-linear equations. The method has ky n applied to the e- 
guation which describes finite amplitude plane wave propaga- 
tion in a viscous fluid in order to compare the results to the 
Beeatctions of a known perturbation solution; the ultimate ob- 
jective being to utilize the technigue for exact solution of 
the corresponding spherical and cylindrical wave problems. 
The first step in this process has been achieved; that is, 
numerical solutions for the plane wave case have been gener- 
ated for a range of values of the various viscous and non- 
linear parameters which are consistent with results obtained 
analytically. Further it is shown that solutions generated 
through the application of parametric differentiation may in 
Peete have greater validity for certain ranges of viscous and 
non-linear parameters. 
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NOTATION 


the coefficient of the first-order term in an assumed 
liquid pressure-density relation 


lee elie Varidiole COcELIClents Of See oe soe 


: : and gd, , respectively 


arbitrary coefficients in a generalized boundary con- 
dition 


small signal attenuation ee = vbk*/2¢ for 
tiqgusds, = vil + (y- 1) /Pr]k? [2C. for gases 


the coefficient of the second-order term in an assumed 
liguid pressure-density relation 


a viscosity number for liquids which represents shear 
viscosity and a phenomonological bulk viscosity 


the parameter of non-linearity = 1 + B/2A for liguids, 
= (y+1)/2 for gases 


etc. - coefficients of the points of the finite-dif- 
ference cube 


small signal sound speed 

specific heat at constant pressure 

specific heat at constant volume 

interval or change in the variable which it precedes 
acoustic Mach number 


Neumann factor 





ie Egerirerentscwor shear and bulk viscosity respectively 

E eae variable inhomogeneous term of the parameter space 
equak lon 

g parameter space variable = d€/dy 

Gr Gt Ives 5,7 ete. - notation for ag/dx, gy os ogy ox cee 


oy ox = ete. 


Ge, a Veo STOrermescevariables for Che 1th value of pp, Ue 


a, 4 manctal NOtation fOr Ehe Cwo-dimensional solution grid 

r an indicator of the strength of the non-linearity rela- 
tive to dissipation = Bek/a 

Y ratio of specific heats = oot 

I. RmOonict nod bec Bessel Tunction Of the first kind 

J, Momeuceil becoel, EUNnCEION Or the first kind 

k acoustic wave number 

K Pekm product 

Bl } a general linear operator 

r wavelength of acoustic signal, also used as an arbitrary 


real number in VonNeumann stability test 


i, 0 indicial notation as in Sn.n’ corresponds to g, : 
, f 


N[ J a general non-linear operator 
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>, ¢ 


Kmenarie viscosity = n/P, 
particle displacement amplitude = U6 74 


particle displacement; subscript indicates the base 
solution, 1=0 


an expression for a calculated particle displacement 
Ole Ope U te 


angular frequency 


constants in the fluid equation of state, for gases 
a — Pot O=FO- ror Jaquids, P= ese, Q=P- Pos 


and y is determined empirically 
instantaneous pressure 

ambient pressure 

Prandtl number 


an arbitrary dependent variable, subscript indicates 
the base solution 


instantaneous density 

ambient density 

condensation = P-P/Po 

a nondimensional distance = fekx 


Shock thickness 
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temporal coordinate, characteristic time 


an arbitrary complex number in VonNeumann stability 
esi 


boundary values of 9 

particle velocity, subscript indicates base solution 
particle velocity amplitude 

discontinuity distance = 1/fek 


a generalized spatial vector, the subscript indicates 
the boundary 


Spatial coordinate, characteristic length 
viscosity number 


parameter of interest = Vb/CoA, subscript indicates 
base value 


u/U. 


a timelike variable in the Burgers equation solution 


a non-physical variable used to transform the Burgers 
equation to a heat equation 





a. INTRODUCTION 

The propagation of acoustic waves in fluids and certain 
crystalline substances has been the subject of long and ex- 
haustive study. Due to the inherent simplifications which a- 
rise in the examination of plane wave propagation, this prob- 
lem has received the closest scrutiny. In the general linear 
theory of sound propagation, the assumption must be made that 
particle velocity amplitudes are infinitesimal, allowing the 
reduction of the problem to the classic wave equation. How- 
ever, in many of today's acoustic systems this assumption is 
no longer valid, and solution of the fully non-linear wave 
propagation problem is required for analysis and prediction 
for practical systems. It has long been recognized that sound 
waves whose particle velocity amplitudes are not infinitesimal 
have phase velocities whose magnitudes change with the local 
Merticle velocity, a result which is not predicted by the lin- 
ear theory. This phenomenom received some attention from the 
philosophers of the nineteenth century, and with the exception 
of the efforts of two researchers of the pre-war decade was 
not treated until the last two decades at which time signifi- 
cant interest was generated because of the increasing sophis- 
tication and power levels of acoustic systems. 

The sOlutions of the 1930's form the basis of this re- 
search. Fubini [1] derived an implicit solution to the plane- 


wave problem in an inviscid fluid and reduced it to an explicit 
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solution for low Mach number. Fay [2], the other researcher 
of the 30's, solved the same equation with viscous effects in- 
cluded. Blackstock [3] noted that though both solutions were 
entirely correct, they were restricted in their spatial re- 
Moms OL applicability. Applying weak-shock theory, he demon- 
strated the manner in which the two solutions are related and 
developed a single function which describes the solution for 
all space in the inviscid case. Keck and Beyer [4] performed 
a perturbation analysis which, for a specific range of para- 
meters, described propagation in the viscous case. Black- 
stock [5], using approximations which he demonstrated [6] to 
be equivalent to those used by Keck and Beyer, was able to re- 
duce the governing differential equation for plane wave propa- 
G@aeron in a thermoviscous fluid to the Burgers equation for a 
boundary value problem, which yielded an analytic solution for 
all space and time. Propagation in relaxing fluids has re- 
ceived moderate attention [7,8], but has generally been a- 
voided. Blackstock [9] has derived analogous inviscid solu- 
tions for cylindrical and spherical waves, but as fate would 
Meme it, an analytic solution for cylindrical and spherical 
waves in a viscous fluid hasnot yet been obtained due to the 
fact that there exists no apparent simple transform, such as 
that utilized in the plane wave case, to obtain the Burgers 
equation. 


There are various parameters which describe the proper- 
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ties of the medium and two which relate to the source excita- 
tion, and taken together they give one an indication of the 
melative strengths of non-linear effects and viscous (dissipa-~ 
tive) effects. The parameters of the medium are vy, the kine- 
Matic viscosity; ¥, the viscosity number; y, the ratio of spe- 
cific heats; Pr, the Prandtl number; 8, the parameter of non- 
linearity; and Cor the small signal sound speed. The expression 
ie a 1) is a collection of thermal and viscous coefficients 
which together with kinematic viscosity relate the thermovis- 
cous nature of gases. In liquids, such an expression is not 
readily determined and is generally replaced by a single term, 
[i B is equal to (y + 1)/2 for perfect gases and to 1 + B/2A 
foeerluids of arbitrary equation of state where A is the coef- 
meerent Of the first order term and B/2 is the coefficient of 
the second order term in an assumed pressure-density relation. 
Analogous parameters for dielectric crystals allow them to be 
treated in a manner similar to that for nieunondecous fluids 
[10]. The two parameters related to the source excitation 
are Use the peak particle velocity, and w, the angular fre- 
quency of the source. Much of the previous effort in this 
field has made use of a parameter TT = Bek/a to express the 
relative strength of non-linearity and viscosity, where e€ = 

5! the acoustic Mach number; k = w/Cos the wave number; 
and a = vbk*/2C,, the small signal attenuation coefficient. 
r~ 1 may be considered as the borderline for the inception of 
important non-linear effects [11]. 


Another important parameter is the discontinuity distance 
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of the particle velocity waveform becomes negatively infinite 
due to non-linear generation of harmonics, that is, X indi- 
cates the point at which a shock first forms. Shocks were 
first predicted by the inviscid theory, and the discontinuity 
distance which has become a reference point in most non-linear 
maeery, was predicted by the Fubini solution, which, in fact, 
is not valid beyond the discontinuity distance. Shocks form 
because the finite amplitude of the particle velocity causes 
the compressive portion of the waveform to 'catch up' to the 
rarefactive portion of the waveform; because a multi-valued 
waveform is a physical impossibility, a shock must form unless 
dissipation reduces the wave amplitude sufficiently. Figures 
la and lb show graphically such a waveform. The utility of [. 
has been extended by Fenlon [12] to indicate whether and where 
a shock will form in the viscous case; for plane waves [ > 4.5 
indicates that shock formation will occur, that is, non-linear 
eetects will be sufficiently strong, so that even with dissi- 
pation, shocks will form. Shock formation marks the end of 
the first of three identifiable zones of propagation of the 
finite amplitude wave; this is a region of strong ee ieinear. 
effects. The second region of propagation is that in which 
men-linear effects initially dominate, but gradually taper off 
as absorption reduces the magnitude of the fundamental and 


even more rapidly the magnitudes of the non-linearly generated 
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harmonics. The limit of this region is reached when the rate 
of decay of the fundamental due to small signal attenuation is 
matched by the rate of decay due to the generation of harmon- 
ics, that is, when the particle velocity amplitudes are once 
again infinitesimal. 

The aforementioned solutions may be discussed with res- 
pect to their relations to the parameters of the problem. The 
Fubini solution predicts shock formation at the discontinuity 
distance, but the solution is no longer convergent beyond the 
discontinuity distance. Fay produced a solution for the most 
stable waveform, a sawtooth, to which the finite amplitude 
wave generates beyond about 3 discontinuity distances. The 
Fay solution contained [f as a parameter, but did not reduce to 
Fubini for [ + ~; Blackstock's bridging function accomplished 
this. Keck and Beyer's perturbation solution is valid for 
very weak non-linear effects (relative to dissipative effects). 
Because of the limited number of terms generated, it is valid 
meee! = O(1) or less. Finally, Blackstock's Burgers equa- 
tion solution is valid for all Tf. Unfortunately the series 
which represents the steady state solution is very slowly con- 
vergent, limiting the usefulness of the solution. Finally, it 
should be noted that an inherent restriction on the size of 
the Mach number is assumed implicitly or explicitly in all 

‘ 
cases, with a maximum value of ¢« = .1l. Such a value would be 


stretching the applicability of certain approximations which 
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have been made in each of the above cases. 

In order to dispense with the aforementioned problems 
with the two viscous solutions, and specifically, to essenti- 
ally eliminate the approximations which were made therein, a- 
nother technique must be used which returns to the basic par- 
tial differential equation and solves it without approximation. 
Parametric differentiation is such a technique. Given a non- 
linear differential equation in which there exists a parameter, 
one need only assume that the solution is continuous in that 
parameter to a limiting value for which a known solution ex- 
ists. This technique has been applied in the past by Rappert 
and Landahl [13] to non-linear flow problems and by Harris [14] 
to aerodynamic sound produced by a rotating cylinder in a vis- 
cous medium. VanDyke [15] has linked parametric differentia- 
tion to work he has performed on the extension of perturbation 
series solutions. The application of parametric differentia- 
tion in this case to a problem whose approximate solution is 
known will (1) shed some additional light on the physics of 
the problem, and (2) demonstrate the applicability and utility 
of a new technique for use in studies of the propagation of 
finite amplitude cylindrical and spherical waves in a viscous 


medium. 
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mie = =LHE EQUATION OF MOTION 

A plane wave propagating isentropically in a homogeneous 
viscous fluid may be completely described by the continuity 
equation, the Navier-Stokes equation for irrotational flow and 


gameequation of state. The appropriate equations are: 


Meecontcinurty op = (1 + 3§/9%) p, CZ ay) 
; Bae ) 4 1 3 2 

2) Navier-Stokes Po Jot = -dp/oax + (an + n')d°E/dx* ot 

(22) 

3) State p= P(p/p,)' - 9 [4] | (2.3) 


where p is the density, & is the particle displacement, p is 
the pressure, n is the shear viscosity coefficient, n' is the 
bulk.viscosity coefficient, and P and Q are constants dependent 
on the fluid considered. The subscripted parameters denote 
rest values of the respective quantities. For ideal gases, 

P = Por Q = 0, and y is the ratio of specific heats, oe 
Por liquids, P = ea, where Cy is the small signal sound 
Speed, Q =P- Por and y must be determined empirically with 
the assumed equation of state above. The equation of state 


for liquids is often written in the form: 


p=p_ + a(2tPa) + By2(PPay* 4 o(®Pa)*)... ea) 
oO Po Po Do 
with — ——ceiMnem@ondensatzon. For small values of s, this 
O 
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equivalent to the equation of state shown (2.3). This can be 


shown by writing 


p= P(1 + 2 fay’ - 9 (2.5) 
Po 
and performing a binomial expansion for eae | < ll, the result 
O 
being 
peepee vp(22loy » Ve!) p(eoho~ , 1. - (2.6) 
Po 2 Po 


Then, with P - Q = Po: A is equivalent to yP, and B to y(y-1)P 

and the ratio B/A = y-l. When one uses this form the implicit 

assumption is made that terms of 02-22) *) ana Gqiecaton Ware loc 
O 


important. For any fluid of low compressibility this is 


clearly valid for small acoustic Mach number. (Note that p/p, 
E ee Oe 
eal + ae and NE e). 


Returning to Equations 2.1-3, p and p may be eliminated 
tO yield a single equation in the particle displacement. From 


Equation 2.3 


op . PY (@)' 30 
OX Po (SO) OX ae 
. , op _ a°F 
Sa@ostituting Equation 2.1 and its x derivative, Ree Dae c 
(l + g5)~? maiicOmbanation 2.7, one obtains, 
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aes 
oP LL . Baar 
ax PYG + Oy) ox? oo) 


Pmoestituting in the Equation 2.2, rearranging terms, and in- 


serting P = 5 OR one obtains the following equation: 
= at oun ae eg Oe 
oo) ED Soe «So Bx? We 
' 
with v = N/Poe b= = + ae At this point one may leave the 


memation as it stands and apply parametric differentiation 
directly or make one additional simplification which resorts 
to the assumption that the displacements though finite are 
still not large enough to require the retention of cubic terms 
in particle displacement. Since previous effort has been dir- 
ected along the second path and since this previous work is 


the point for comparison to determine the validity of the ap- 


proach, this additional simplification is made. (1 + ge) 7 y+ 
is expanded binomially to yield: 
DE. = oe Gan (28) 3 
coat =) =e Caml ema ypele ae oa) * + ol 
C20) 


The first two terms only are retained since this factor multi- 


plies the right hand side. The resulting equation is: 


Gta Vor ek Yous _ 9&6 907E 
ax? ~ C7 at? + COZ oxzae ~ '*8) GN OxF Sher 


O 
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where 6 = (y+1)/2 for gases and § = 1 + B/2A for liquids. This 
is in essence the equation treated by Keck and Beyer [4], a- 
feng with the boundary condition £€(0,t) = Se QcOsuE, in. Enea: 
perturbation analysis, after terms of O(a/k) 4) are eliminated. 
The Keck-Beyer solution is computationally simple, but it is 
severely limited by the labor required to calculate the terms 
Of the perturbation series. In fact, according to the auth- 
ors, the perturbation series 1s probably only convergent for 


Mee (i-e 29 


)f/4a < 1, which restriction corresponds, for exam- 
ple, to an acoustic pressure of .125 atmospheres at 1 MHz in 
water; one could easily argue that this pressure hardly qual- 
ifies as finite amplitude,but then if one is able to measure 
finely enough, any acoustic signal can be considered to be of 
finite amplitude. Actually, the series has been demonstrated 
to be everywhere convergent but is still limited in applica- 
bility, by truncation, to values of the parameters as defined 
above [7]. It is to be noted also that fek/4a = [/4; thus the 
Keck-Beyer solution will not allow [ > 4, and since [ > 4.5 is 
required for shock formation in the viscous plane wave case, 
this solution deals with non-shocked waves only. 

Blackstock has Ferenc ced that the governing equation 
may be reduced to a Burgers equation by applying the same 
assumption utilized in the perturbation analysis. This as- 
sumption, in fact, is the same one which allows one to arrive 


at an analytic steady state solution of the following problem: 
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evo e. 9.2 97E 
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namely the small signal problem for plane wave propagation in 
meeerscous fluid. This assumption, that a/k << 1, is valid for 
most fluids, with the exception of the most viscous ones, at 
mmyeceasOnable acoustic frequency. In any case, the effect of 


making this assumption is to allow one to utilize the expres- 


Sion: 

Bees | of 

ax ~ Got’ ee 
mereaducing the full equation to a Burgers equation. This 


expression is, of course, the differential equation which des- 
cribes the propagation of a small amplitude wave ina lossless 
medium. The reduction procedure follows. Starting with Equa- 


tion 2.11 (an approximate one which has been obtained by as- 


Suming a small Mach number), one first substitutes ee. = Ey 
O 
Some the right hand side to obtain: 
Ral vb _ _ 28 
Sx on Ste t Che Sext ac. Stoxx eae’ 


Then note the following identities from Equation 2.13: 


eee 7 
ry Ce SEt (a) Ca eee Ke 
ao. 2 ae = 
See eet 7 et (b) 6 See 7 Souk So eee (d) 
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= th 2 
oe a oe St (e) ae CF Set (ce) (2.15a-f) 


Substituting (d) into the viscous term of Equation 2.14, 


vb _ vb Eee epee 

eo ee Tom ee ee em (epee 7 “ofxee? and 
re Vb _ Eee 

oxx oe Stet 7 OC? Sete 7 Sobxtt? = Co ep otea aa 


Substituting (e) and (f) into the non-linear term of Equation 


eel , 


NO 
W 


B 


7 e- Se oxy = = Wee 7 eos 
48 i; a 
mC. (a eee G2 Ee cee? 

O O O 


The equation after substituting and manipulating signs be- 


comes 


Ht vb 


2 VV = ie uy = = 
OEE Z OES tt ee Coe xx oe ek one! 


(72 aly 


Performing an operation which is the integral equivalent of 


differentiating with respect to the operator a - Coat, one 
Obtains 
2 " vb 3 Po 2 
Coe t 2 See F Coby = 7 BOOS coeee 





ei 


Differentiating with respect to 't' and substituting u = Gy 


2 vb 3 = 
Cc Ge ae BC ou Up (2.19) 


Cet Set ox 


Finally, with the variable transform 


Se — iE — te a = t 
Oo 
one obtains 
il 
Cae = BC ou Ur =F 5 VDup rg ’ (2220) 


which is Burgers equation for a boundary value problem. The 
Significant point here is that one is considering an equation 
which is non-linear, and which includes viscous terms. The 
approximation a/k << 1, is applied two times to reduce the 
equation to Burgers form, each time to simplify viscous and 
non-linear terms. This is indeed an acceptable engineering 
approximation; however, it certainly is desirable to have the 
ability to solve the equation without resorting to this ap- 
PEeoximation, particularly for highly viscous fluids and per- 
haps certain dielectric crystals with large stress coeffi- 
cients. 

As noted by Blackstock [6], the Burgers equation solu- 


tion has one further limitation which might be considered more 





a 
severe. The steady state solution to the Burgers equation 
1s 


aoa 


ee 
v = hs En! i) ee a) cos ny 1D Da 


where t is related tou by, 


_ 2 
a r [log o] (22 2)) 
Here, oa is the Neumann factor (E = 1, all other A 2) IT 
is the nth order Bessel of imaginary argument; o = fekx, a 
non-dimensional space variable; and y = wt', a non-dimensional 
temporal variable. For values of [ > 50, the series represen- 


tation of the solution is very slowly convergent, and there- 
fore difficult to use practically in numerical analysis. Con- 
mementiy, £or ~ > 50, Blackstock utilized an asymptotic ap- 
meeach to the solution, which in fact is not valid near the 
mergin. TO solve the Burgers equation for Bekx << [, he re- 
sorts to a perturbation analysis in nae but recalling T= Bek/a, 
this is in some sense similar to the Keck-Beyer approach which 
utilized € as the perturbation quantity to solve the same basic 
equation. Blackstock's solution for small x contains terms of 
first order in the perturbation parameter; for I >> o, this 
should indeed be appropriate, however the coefficients derived 


contain an infinite series of Bessel functions which may be 
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slowly convergent themselves. Herein lies another motivation 
for the parametric differentiation solution, the desire to 
produce a unified solution applying the same assumptions and 
having the same limitations throughout. 

In any event, there iS no simple transform which will re- 
duce either (1) the full equation valid for high Mach number, 
or (2) cylindrical and spherical finite amplitude wave equa- 
tions to the Burgers form. Therefore, an alternative ap- 


proach to the problem is required. 
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fol. PARAMETRIC DIFFERENTIATION 

Parametric differentiation iS a procedure which allows a 
non-linear equation which involves a parameter to be solved by 
transforming the equation to a linear equation in parameter 
Space. This is particularly valuable for cases in which nu- 
merical solution will be required because it eliminates the 
possibility of obtaining multiple roots in solving non-linear 
difference equations. It is essential to note at this point 
that some insight into the nature of the solution at the ex- 
tended range of the parameter of interest would be extremely 
helpful. If the solution is singular in the parameter, that 
is, if the nature of the solution changes radically as one ap- 
proaches the limiting value of the parameter in the known base 
solution, then the technigque is not applicable. Thus, para-~ 
metric differentiation would never be applicable to the entire 
class of singular perturbation problems. The method of appli- 
@2ae10n follows. 


Suppose one is given the general problem: 


Nig (X, ey WV) ] = 0 (3.2) ) 
Ol(k,, ts w)) = 6, (t; wp) (3. 1a) 
1 


OL (Xp e ee coe (3b) 


ct 


abl (Kpe t; VI] + co (Ry, ts P) = 0, (try) (3.1c) 
3 





where N is a non-linear operator; X is a spacial vector coor- 
dinate; t is a time coordinate; is the parameter of inter- 
See, and the 0's are functional values of 9, ee or ag + ch. 


at X_ as appropriate to the problem, n indicating the normal 


B 
derivative. Initial values could also be prescribed if appro- 
emetate. If a solution 6 = Geet vo) exists at some limiting 


value of the parameter Wor and it is desired to obtain a solu- 
tion to the same problem with identical boundary and/or ini- 
tial values at Be + Ay, then parametric differentiation may be 
applied. First, differentiate the governing equation and 
boundary and/or initial conditions with respect to the para- 


meter to yield: 


mig (x, t; py] = 0 (3.2) 
g(X,, ti b) = ha y) (3. 2a) 
g,(%yr ti v) = 9g, (ti v) (3. 2b) 
2 
ag(X,, ti ¥) + cg (X,, ti ¥) = Tp (thw) (rc) 
where 
ee a = 
etx, ee py) = av (> (X, ey W)) (353) 


One should realize that this procedure cannot be applied 1f 
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the parameter appears in the basic equation raised to a power 
less than one, and if the limiting value of the parameter Yo 
= O. 

The operator, L, is a new linear operator obtained through 
Gifferentiation. The new equation with associated boundary 
and/or initial conditions may now be solved in parameter space, 
either analytically or numerically, for successive values of 
the parameter ~. As the equation is solved for each value jy, 

a quadrature is performed which solves Equation 3.3 and re- 
turns a value or functional form which, when applied to the 


solution for Vs yields the desired solution at Ws. There 


-|’ 
1s no numerical restriction on the range of wW; however, there 
certainly may be a physical limit, and some insight is re- 
@mared here to avoid blind application of the method to obtain, 
for instance, solutions for a range of wy which is beyond the 
region of validity of the basic equation. The quadrature may 
be expressed as follows: 

ime = G4 1X, t7 v.45) + | sav (3.4) 
For more complicated equations which require numerical solu- 
tion, this step can be quite complex, depending on the form of 
the g versus y~ curve. In some cases, trapezoidal integration 
may be sufficient; in others, higher order integration techni- 


ques may be required. It is clear in examining the form of 
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the equation above, that this is, in essence, a recursion type 
relation and consequently it demands a known initializing 
mumeti1on for i=l. 

ie chis work, Ehe initializing function is the Fubini 


solution to the problen, 


3g 28 g7E _ 2 ae 
mY ox’ ste ~ “o ox? Sa 
Men the boundary condition £(0, t) = - EZ cosut, this being a 


duplicate of Equation 2.9, without the dissipative term. Thus, 
the parameter of interest here will be some non-dimensional 
guantity which contains v, and the quadrature will be per- 
formed with respect to the non-dimensional quantity. An ex- 
plicit form of the Fubini solution was given by Keck and 


Beyer [4] in terms of particle velocity, 


ee) J (nkx) 
ect ait. I ee eit =i) (3.6) 
n=L 
where em Pome heownetmOracr Bessel fLUunCcELOnN and Kk = Bek. This 


Semicion is limited in Mach number (AO <u) eSinee. ao binomial 
Oo 


expansion in Mach number was utilized to arrive at this result. 
An exact implicit solution is available which could be used for 
analysis at greater Mach numbers, 

= _ Wx Slee 
u(x, t) = Uosin (wt (Tee a 
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However, the numerical problems associated with an iterative 
solution at many points would be difficult to resolve. This 
solution predicts the formation of the shock at the point 
where the velocity profile becomes negatively infinite, that 
is, where ou/ox > -~, Differentiating Equation 3.7 with res- 


pect tO xX, 


oy OE - (5) 
dU _ WX ie, eS ere ye 
\e U,cos (wt G (1 + 5 a) Dl a (1 + 5 mel 
Oo Oo O Oo 
ae ee 
vy+l wx Ve Ul v=o Veo 
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ce 
O 
This quantity becomes infinite when 
2Y 
zs yori wyy-l ytli Yq w 
x (1 + ae ) / EEG cos( ) (3210) 
Oo on 10 
Noting that 6 = ee Gon ee k, and that the first shock 
Oo O 


will occur at the head of the waveform where cos( ) = 1, and 





Oe 


meee that uw = 0 at the shock, the discontinuity distance X 
i 


ma ek iemOo-aine dua sicethe lamit Of convergence of the ser- 
Hes SOlution stated above. Thus, using the Fubini explicit 
solution as a starting solution, one is limited to application 
of the technique of parametric differentiation in the region 
O< x < X and for Mach numbers which are small compared to 
unity. It is clear from this discussion that one is fundamen- 
tally limited when applying parametric differentiation by the 
mality Of the base solution. It is unfortunate that in this 
case, the base solution is limited with respect to Mach number; 
on the other hand, signals of sufficiently large amplitude may 


no longer be thought of as acoustic waves and may propagate in 


a different manner. 
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ive APPLICATION OF THE METHOD 

Tt is one thing to discuss the use of parametric differ- 
entiation in abstract terms, and quite another to execute its 
aeplication, there being no real ‘theory’ with respect to how 
memshOould be applied. In any case, the first step in the sol- 
M@tion technique is similar to that for most other numerical 
solution procedures, that is, non-dimensionalize the equation. 


Starting with the reduced equation 


il Vb = 
bx Cae eon Cees = BE ex (4.1) 
one introduces characteristic lengths and times. This being a 
wave propagation problem,it is reasonable to select aaa X and 
t. = = , where dX is the wavelength; then the non-dimensional 
O 


Variables are: 
x = x/X = x/h t= t/t. = tc /r and € = E/% = E) 


Then, Equation 4.1, in terms of dimensionless variables, is: 


523 tt * CoA °xxt - ae OR oo 
This equation is as one would expect for a ‘'good' non-dimen- 


Sionalization; the coefficients of the first two terms, which 


in fact should dominate, are of O(1). Alternatively, with the 


=o 


emgation in the form 


) = C7E (Aaa) 


28 
(1 + Ex) Sone . eee Oe 


a similar non-dimensionalization will yield 


Dee 28 ee = Vb ee _— PA 
(1 + 53) (Ef CoA ee ~ oe oo 


which demonstrates that the relatively large numerical value 
O£ the coefficient of the non-linear terms in 4.2 is simply a re- 
sult of the expansion performed, and should not be considered im- 


proper. The viscosity coefficient appears as a which 1s 
O 


like an inverse Reynold's number. However, because Cx 1s a 
propagation speed and not a particle velocity, this 1s not a 
classic acoustic Reynold's number. In any case, this quantity, 
which will be termed 'wW', is the parameter of interest, on 
which the quadrature will be performed. The solution for ,p=0O 


1s known; it is the Fubini solution. Then one has 


u(x, t: yp.) = 20. J Un) cinn(wt - kx) (4.5a) 
O O a Nnkx 


or, in terms of particle displacement 


_ _5Ug & Jy (nkx) - 
eek, tr ye) = aD Ly cosn (wt - kx) (4.5b) 


Wwetrech corresponds to the function 6 imenewtReoret ical devel 
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opment of parametric differentiation, and with =. = ae 
O 
The next step in the development is to differentiate E- 


guation 4.2 with respect to yj, yielding 


Ween Ieee tf = 2RE Og, + 2 RE LoL (4.6) 


where 


Q 
TY 


(457) 


Q 
lI 
= 


and the umlaut notation has been dropped from the non-dimen- 
Sional variables. One would immediately look at this equation 
and express horror at what has been done; there seems to be 
added complexity rather than simplification. However, one can 
now treat the equation as an inhomogeneous, linear equation 
Watch variable coefficients, which should make it more tract- 
able for numerical solution. The only non-linearity which re- 
Mmeens 1S found in Equation 4.7. It should be noted that the 
non-dimensionalization in terms of the size of the coefficients 
has not been adversely affected. In fact, the dominance of the 
first two terms has been enhanced since by! cot S and e« << l. 
Alternatively, differentiating Equation 4.4 one obtains 


26-1 28 
ete) eG (ene > Ye.) t (1 +e) tg.) 
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xxt )=g (4.8) 
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which really looks as if this "simplification" could have been 
done without. However, if one multiplies through by (1 + oe 


the following equation is obtained 


(28) (1 +E) Pe, - ve ag + 1 + BPM ig - vo, 
Sxxt? as! bx) Sax oe 
2BE Gg + (1 + bP ig - We - Se) = (bt ED, 
(4.9) 
which is of the same form as Equation 4.6. This now is the 


equation which will be utilized to solve the finite amplitude 
plane wave propagation problem without approximation, save 
those which may be necessary to obtain a base solution. In 
this case, this allows a significant extension because this 
equation does not require the approximation a/k << 1, and 
therefore, the Fubini solution may be used as a base solution 
for the problem of propagation in highly viscous fluids since 
it is limited by Mach number only. 

The quadrature required to extract the desired solution 
is essentially dependent on the problem itself. Landahl and 
Ruppert [13] employed Runge-Kutta integration to solve a boun- 
dary-layer problem. In many cases, it is not necessary to re- 


sort to such exotic techniques; in this analysis, g varies quite 
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slowly for small values of jp; a rapid non-linear change in g 
with WwW only occurs as becomes very large, that is, as the 
viscosity becomes high. PoEMnOstetluvas, Including) yery 
viscous ones, a value of w sufficient to produce rapid varia- 
tion of g with yw would not be attained at any reasonable acou- 
Stic frequency. Consequently, the integration technique em- 
ployed in this analysis was a very simple straight line inte- 
gration, which required one iterative step. The initial sol- 
ution of the equation for g on the ith step was computed based 
on the ith value of wW and a solution Ei (x, t) which was ob- 


tained in the following manner 


(4.10) 


Ga ee) = B53 bx tt) eat G;_1 (x, ier (ps = v5 4) 


1 
With this new solution, gi (x, Ea Modirlcatron OF Ei (x, te) 


was performed as follows: 


iu 
yO t) = Ey bey t) + Z(Gy (xr ED - oy x, t)) 


E (Ys a W541) 41 1} 


The net result being: 
E.(x, t) = €,_,(x, t) + $(g, (x, t) + 9,_, (x, ¢)) 
age ieee 2! Ze ih” 


v= Vas) (4.12) 





which is in essence a trapezoidal integration. Clearly, this 
scheme would not be useful for g functions which vary rapidly 
and non-linearly with yw, that is, if Equation 4.7 exhibits 
highly non-linear behavior. This scheme could be used in an 
iterative manner, however, with continuing refinement of the g 
values at each step. The entire process used to arrive at the 
extended solution for &(x,t) is shown schematically in Figure 
2. it may be noted that there is a requirement to obtain var- 
iable coefficients at each step, this of course is a situation 
which is not enviable in a numerical solution in which the co- 
efficients cannot be described analytically, particularly when 
these variable coefficients involve derivatives of the depen- 
dent variable. In this case, all the coefficients involve de- 
rivatives. A simple solution to this problem has not been 
found; however, in this analysis, numerical differentiation of 
the € solution to obtain the coefficients does not seem to 
have affected the results. 

Another not so simple problem in the application of para- 
metric differentiation is the treatment of the boundary condi- 
tions. One must be extremely careful in not trying to read 
meemmuch into the Te In this analysis, two approaches 
were considered for determination of boundary values for appli- 
cation in the finite difference scheme, the result being that 
the most natural and straightforward approach was eventually 


ehosen as the proper one. In analyzing this problem, one can 
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make some subjective evaluations of the expected behavior of 
Bepoeased On what the gAy product represents. Dissipation will 
tend to slow the non-linear generation of harmonics as the 
wave propagates, and if dissipation is strong enough, one 
would expect the waveform to remain sinusoidal, if such is the 
nature of source excitation, no matter how strong the non-lin- 
earity is in absolute terms, that is, for [ < 1 the waveform 
Smemlad not steepen appreciably. In any case, near the bound- 
ary the effect of gAw should be to decrease the amplitude of 
the signal. Aw is always positive, therefore to obtain such 

a reduction in amplitude, one would expect the sign of g to be 
opposite that of the base solution. Attempting to produce an 
analytical model which fulfills one's expectations near the 
boundary may be self defeating. In this analysis it was ini- 
tially considered desirable to introduce a small perturbation 
of the boundary condition in terms of ww, and consequently the 


boundary condition was rewritten as, 


{1} 


moO, t) = = +0) cos (wt) (413) 
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which yielded the intuitively desirable result in terms of 


SGN, 


3 Es 
5p 6b (0, Pe (O78) = sae) zI cos (wt) (4.14) 


This model also seemed appropriate since it produced a bound- 
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ary value of g which decreased slightly as w increased, corres- 
menaing tO an anticipated increase in the viscous effects. If 
One recalls the nature of the quadrature, the only effect of 
this model would be to force the solution near the boundary to 
mice, on the desired form. This ‘forcing’ of the solution by 
injecting one's subjective expectations into the modelling of 
the boundary conditions was not particularly bad, but neither 


micett accomplish anything with respect to improving the solu- 


tion. If one instead blindly applies the ah Sperator Luo rane 
meunadary condition, the result is 
3 = 
oO, t) = BY ~oCOs HE) ) =0 (a2s5)) 


Whach does absolutely nothing to convince one that the proper 
eeeeeOor Quadrature will occur near the boundary, but instead re- 
lies on the differential equation to properly predict the g 
values everywhere. In the final analysis, this was the cor- 
rect procedure; there were in fact very small differences in 
the solution values obtained with the forced boundary condi- 
trons, and those obtained with the ‘'natural' condition. 

One may wonder why initial conditions are not important 
in this problem since it is essentially hyperbolic in char- 
acter, and such problems are generally thought of as initial 
value problems. The answer is that when one is interested in 


the steady state solution only, for a problem with a periodi- 
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Seably varying boundary condition, and in this case that is the 
Petition desired for &(x, t), one may set the condition that 
the initial values occurred at some sufficiently large time in 
the past that all start-up transients have disappeared, and 
therefore only the boundary conditions determine the behavior 
eeene solution {16]. It is only natural to expect this physical 
interpretation to be carried over into the parameter space ap- 
plication, and consequently, initial values have been of no 
concern. 

Tt is valuable to note that the final result for the cho- 
sen maximum value of p need not be the only output. For a 
maven Be product, a solution for any IT consistent with the ac- 
curacy obtainable (as [ + ~, the accuracy required to obtain 
meaningful results is increased for a given Be because the 
differences between the inviscid and viscous solutions go to 
zero), may be obtained by selecting the appropriate yp = pe and 
simply calling for output at the end of the second iterative 
@eeee FOr a given fluid this would allow one to analyze the 
losses at different frequencies for a given fe and, in this 
Stee, Out to the discontinuity distance. In fact, tabulations 
of g versus wy) at specific points could be stored so that the 
solution could be extracted in the future simply by integra- 
mmegeand applying the correction at the appropriate point in 


Eme base solution. 


It is useful to compare the base solution waveform with 
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the anticipated extended solution waveform, the comparison in 
this case being very simple because the extended solution or 
a very good approximation thereto iS well known. As it de- 
velops, the incremental differences between the base solution 
and the extended solution are expected to be in phase with one 


a@mother, if they are written in the form 
Helx, t) = Ee lx, te) ee EG xs t) (4.16) 


where A&(x, t) 1s the difference between the calculated solu- 
tion and the base solution when integrating Equation 4.7 from 
Ye Eo Vee The observance of A&(x, t) varying in the expected 
manner is a key clue to the success of the application of the 
method. This behavior was in fact observed in this analysis, 
and will be discussed in a following section. Once one has 
feeeloped an equation in parameter space and made a subjective 
evaluation of what he expects to be the functional form of g 
and of A&é(x, t), he is ready to proceed to the numerical tech- 
nique to be used for solution. Ideally, it would be hoped 
that the transform to parameter space would lead to an anal- 
fee Solution of Equations 4.6 or 4.9 and Equation 4.7, or at 
m@emvery least, Equation 4.6 or 4.9. Such is not the case 
here, and numerical solution of both was required. Thus, to 
analyze the results, and prove the worth of the method, a sub- 


jective feel for the results was required. 
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V. FINITE DIFFERENCE EQUATIONS 

After the equation of motion has been derived, non-dimen- 
Sionalized and its analog in parameter space has been obtained 
Pmmenugh dirrerentiation, it may still be a significant task to 
obtain the solution of the linear g equation. In this case, 
Pamation 4.6 or 4.9 describes the function g with boundary 
Pemaition g(O, t) = 0. In analyzing the equations, it is noted 
that they are essentially wave-like in nature, at least close 
to the boundary, since the Sey and Sig terms dominate. Adopt- 


ing the approximation ae —&, (in non-dimensional form) and 


c 


differentiating with respect to x, one notes that oer atl 
and since nee 2° at the shock in the inviscid case, then oe 
>+ eo at the shock. This wave-like behavior may not then be so 
readily observable away from the origin since the coefficients 
involving oe may become very large in the regions in which 
the shock is forming, or has formed in an extended solution, 
in which case the g, term becomes important. The or term 
remains small since wy is generally much less than one. The Ge 
term, of course, is the one which results from the inclusion 
of non-linear terms, and it is entirely consistent that it be- 
comes important in the regions of the waveform which exhibit 
the result of non-linear propagation, namely in the region at 
the head of each compressive portion of the waveform and at 


me tail of the rarefactive portion. Once a shock has formed, 


it is conceivable that the J, term will alter the character of 
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the equation to the extent that it is no longer wavelike; this 
is not of immediate concern since this analysis is restricted 
to the first region of propagation. Therefore, in attempting 
to solve this equation, it was treated as if it were wavelike 
throughout. As previously noted, the literature is rich with 
various methods of solution of the fundamental equation of mo- 
tion, techniques which, if applied judiciously and with essen- 
tially the same approximations as those utilized to solve the 
equation of motion, could possibly be applied to Equation 4.6 
or 4.9 to yield analytic solutions. However, because of the 
form of the variable coefficients, this analytic solution 
would at best be extremely difficult and most likely would be 
impossible; so in this analysis solution, by finite differences 
of the parameter space equation was accomplished, and numeri- 
Gal integration of Equation 4.7 was performed. Since the e- 
Gimeeaon 1S essentially hyperbolic in nature, a finite differ- 
ence scheme appropriate to a hyperbolic system was the obvious 
enoice, 

There are two broad categories of finite difference 
schemes for partial differential equations, which are termed 
Smemicit and implicit. Adopting a standard notation for a 
nine point difference mesh, Figure 3, an explicit scheme, is 
One which has only one unknown in the i+l column or the j+1 
row, depending on whether one is "marching" in the 1 or j dir- 


ection. An implicit scheme is one for which there are two or 
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more unknowns in the itl column or j+l row. To solve the ex- 
Piepcit scheme for a single value in the future (marching dir- 
ection) one need know only the values of the dependent varia- 
bles at all present and past positions on the mesh. The terms 
future, present and past apply in a spatial as well as tempor- 
al sense. Thus, for the explicit scheme, if one has initial 
Gonditions valid for all space, or if one has initial condi- 
tions for a half space and a boundary condition, a complete 
solution can theoretically be generated numerically for a pro- 
perly-posed problem. The implicit scheme requires solution at 
two or more points based on present and past points. Thus, 
one generates two unknowns in the first application of the 
mesh, and an additional unknown with each subsequent applica- 
tion, yielding N+l unknowns with N equations. Thus, an addi- 
mremal condition iS required in order to apply an implicit 
scheme; generally it is a boundary condition at some point in 
Space. In any case, this requirement essentially restricts 
the use of implicit schemes in infinite domain problems such 
as this one. 

Having been limited to the use of explicit solution 
schemes, there exists one additional problem of significant 
mipert; that is, stability of the scheme. Unfortunately, it 
memattficult to consider stability until one has actually de- 
veloped the finite difference equation, stability of the scheme 


being a function of the way the scheme itself is developed. 
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Proceeding to finite difference approximations of the Equation 
4.6 or 4.9, one finds a literature which is rich in references 
to the wave equation for infinitesimal amplitude waves, in a 
lossless medium, but which has very little in the way of ref- 
erence to non-linear equations or equations for a viscous 
medium. Recalling that the parameter space equation is essen- 
tially wavelike as long as the coefficient of the d,, term is 
relatively small, one can use the finite difference approxima- 
tions to the wave equation as a Starting point. Since an es- 
sential aim of this study was to demonstrate the utility of 
the method of parametric differentiation in non-linear wave 
problems, a second order finite difference scheme was selected; 
it is clear that greater accuracy could have been obtained with 
a higher order scheme. The following finite difference approx- 
imations were used for each of the derivatives of g in either 


Equation 4.6 or 4.9: 


Z 7 ; . 

Gee = (Ggu1, 4 — 294 5 + Syy7 5) (Ax)? + OC (Ax)*) (Sn) 
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where Ax and At denote the step sizes in space and time res- 
pectively. It must be stressed that these approximations, par- 
ticularly Equations 5.3 and 5.4, are not the only ones which 
could be used. For instance, a forward or backward difference 
could have been used in place of the central difference of E- 
quation 5.3. This would, however, have changed the order of 
the approximations. The form of Equation 5.4 is that of a 
second central difference in space and a backward difference 
in time, the backward difference being necessary in order to 
preserve the explicit nature of the scheme. A central differ- 
ence or forward difference in time for Equation 5.4 would have 
not only destroyed the explicit nature of the scheme but would 
also have yielded an unknown to be determined at the point 
ay, j+l whose coefficient, wW, would have been much smaller 
than those of the other terms, and could have led to an in- 
stability, even if one were uSing an implicit scheme. This 
problem could be resolved in an implicit scheme by substitu- 
ting for Equation 5.1 a weighted time average of second order 
approximations of Tyx° 

Representing the variable coefficients of Seen Io 


and g by A Ay, Az, and A, respectively, Equation 4.6 or 


ie 4 


4.9 becomes, upon substituting the finite difference approxi- 


ox 


Mations 5.1-4: 
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aes = F C35) 
where F is the inhomogeneous term. At this point, one may 
look at Equation 5.5 and note that it has one unknown in the 
freerOw 1f One iS marching in space, and two in the i+l column 
if marching in time; so that it is well suited to treatment of 
an initial value problem explicitly, but is implicit if used 
for a boundary value problem. In order to apply the scheme in 
Gamexplicit manner to the boundary value problem, one must make 
Pmeaaitional approximation. To solve for g(it+l,j), one needs 
information in the i-l and i columns, plus information at the 


Point g. Because the third order derivative is present, 


ele) 
there is no way to avoid this problem. However, considering 

the fact that the coefficient of the third order derivative is 
small compared to the coefficients of the other terms, it would 


not be unreasonable to presume that if one had an approximate 


wee Of g. 


ieee Picimanrcolutdeon Of reasonable accuracy could 
f 


be Obtained for g. 


i+l,3! Whewmwould besbased On five pieces of 


accurate information and a single piece of approximate infor- 
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mation. Considering Figure 4, which depicts the finite-dif- 
ference grid with the coefficients of the various points in 
the solution cube annotated thereon, one notes that the ap- 


proximation required at g need be made only once in 


i+l,j-l 


each spatial step. After solving for g. at the 'base' of 


itl,j 


feoangle spatial column, g 1s now known for all subse- 


itl,j-l 
quent applications of the cube as one marches out in time. The 
approximation used to arrive at this "corner point" is a Tay- 
lor's series expansion in which the derivatives at the point 
about which the expansion is taken are approximated by back- 
ward differences accurate to O(Ax). Again greater accuracy is 
obtainable, but was not required because of the relative small- 
ness of the coefficient. The resulting expression for the 
Pee@mner point" is, 


4 Coron 
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It must be noted that this approximation is a potentially weak 
link in the method, due to the fact that near the head of the 
waveform the x derivatives become large, and a truncated Tay- 
lor's series may not be of sufficient accuracy. In any case, 
it did not appear to significantly affect the nature of the 
results in this analysis. 

In the discussion of explicit finite difference schemes 


one must inevitably be concerned with the stability of the 





scheme as it relates to step size. According to the Courant- 
Meredraichs-Lewy criterion for stability of finite difference 
schemes [17], the ratio of the step sizes must be such that 
One does not attempt to find the solution at a point which is 
outside the region of influence of the known data, the region 
of influence being determined by the characteristics emanating 
from the endpoints of the known data. Unfortunately, this 
problem is difficult to address here because of the nature of 
Paemgoverning differential equation. Typically, third-order 
partial differential equations are not treated as such, but 
are instead reduced to some more tractable form; in this case 
the approach using the Burgers equation was such a reduction. 
Consequently, little appears in the literature with reference 
memotability for third-order equations. An unsuccessful in- 
Mesergation was conducted to attempt to arrive at the exact 
Siamacteristics by reduction to canonical form. However, this 
equation is wave-like, and one could reasonably expect the 
characteristics to be similar to those of the wave equation, 
maemeemat the stability criterion would be similar in form. 
Working with this foreknowledge and applying a modified Von 
Meemann [17] stability test, an exact stability criterion will 
be established. The VonNeumann test essentialiy consists of 
examining possible exponential solutions to a finite differ- 
ence scheme to determine when they may grow without bound 


given finite initial or boundary conditions. With this know- 
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ledge one may adjust the stepsize so that the scheme will be 
Stable. Unfortunately, the VonNeumann test is applicable to 
finite difference equations with constant coefficients; in 
this case the coefficients are variable, and the only recourse 
is to treat the coefficients as if they were constant, derive 
Piemstability criterion, and then discuss the implications of 
their variability with respect to how it affects the stability 
eracerion. 

The exponential solutions to the scheme which may grow 


without bound will have the form, 
(Se7) 


where the subscript notation m,n corresponds to the notation 
fmeeiccad previously. i’ is any real number and @ is an arbi- 


trary complex number. Thus, boundary data may be written, 
=e (5ees8) 


and the right hand side may be considered as a typical term in 
a Fourier expansion of the boundary data. Substituting Equa- 


meomeo./ into 5.5, ignoring the inhomogeneous term since it 


eeenot affect stability, and dividing through by ee 


a 


emer obotains 
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+ (-A,/ (At)? + 2A,/ (Ax) 2At)e*” + (-2A,/(Ax)? + 2A,/ (At)? 


2A,/ (Ax) *At) - Ane” / (At)? - age Pe 147 (ax) 2at 


+ (A,/(Ax)? ~ A,/2Ax + A,/ (Ax) ht) e*® = (5.9) 


Rearranging terms such that they are grouped by coefficients, 


the following results, 


10.7 -i6 1A, 7 IK ot 8_71 
A (ae) FR » (ea) — Az ( ie ) 
16. -i86 ~ix 
+e -2, ,i-~e = 
# Rg (Tae) RE) = 0 ed 


Recognizing the various terms as trigonometric forms, one ob- 


eains 


: (4A, sin?3) / (Ax)? + (4A,sin?4) / (At)? - iA,sin0/dx 


: A, (4sin?2) (1-e"**)/(Ax) 7At = 0 fee 1) 


Now solving for the imaginary portion, 


-A,sin@/Ax - 4a,sin?5 SPV Aticore AL = 0 ilo 


which identity indicates that the product, 
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AxAt = -4a,sin?3 sink/A,sin0 (So) 
Sxoula hold, but this is clearly absurd since it indicates 
that as Ax > O, At may become arbitrarily large as well as 
negative. Reason and foreknowledge of the stability criterion 
for the wave equation dictate that this result has no meaning. 
Returning to Equation 5.11 and equating the real parts, one 


obtains 
Z een: 2 sence ae - 20 7 2 
( A,sin 5) / (Ax) + (A,sin x) At (A, sin x) (1 COS) )7 (Ax) ht 
= O (5a) 


Or 


A ieee 
Pn! = é 2 emeie) 
D A, (At) 7+A At ( ) : 


4 T—COs) 
Now one argues that to prevent the assumed exponential solu- 
tion from increasing without bound, @ must have no negative 
imaginary part. However, because of the quadratic form, if 86 
has imaginary roots, they will appear in conjugate pairs; 
therefore 6 must be real. Since A was assumed arbitrary, this 
requires the following restriction, 


A, (Ax) ? 
— A, (At) 74A, At (1-cosh) se (5.16) 
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Now, noting that the maximum value of the expression above is 
obtained when l-cosi = O, it is clear that a sufficient ex- 
pression for the above is, 


A, (Ax) * 
oS 
A, (At) 


[A 
- 


(5217) 


which is of precisely the same form as the stability criterion 


for the wave equation, in which case A, = eee) eta Ss 


case Ay and A, are variable with Ay = je bee and A, a (l+ €.) 


mor Bouation 4.9. Since by 1S approximately equal to the lo- 


2d 


cal acoustic Mach number, which is taken to have a value of 
po -) Peis wanalvsis, the criterion is only slightly dif- 
ferent from that of the wave equation. However, it iS impor- 
tant to realize that the ratio of step sizes cannot be unity 
as one may have anticipated. It is also comforting to note 


that A the coefficient of the oe term, did not enter into 


3! 
the criterion in any way, lower order terms generally having 
little effect on stability. This statement, however, deserves 
further comment. Recalling the discussion of the size of the 
respective terms in the parameter Space equation, the gy, term 
has a coefficient, See which will become very large in the 
region of a shock. It is inconceivable that this term will 
not affect the stability of the scheme in such a case. The 


point is that when bux becomes very large the entire nature of 


the equation does change with the gy, term being the critical 





S50 


term. In this instance the stability analysis is no longer 
valid, and one would be forced to rely on some other stability 
test. This problem is not great here because the analysis is 
mestrictea to the first zone of propagation. Also, the dis- 
cretization of the displacement data limits the maximum value 
Eenat a can attain in the numerical scheme. This is, as yet, 
an unresolved numerical problem; to allow ane GCOvdttain Ls 
true value in a nearly inviscid fluid, an exceedingly small 
step size in space would be required. However, if one consid- 
ers the actual size as IMay attain in a typical propagating 
wave, this problem may be put in a proper light. Blackstock 


[5] has defined a shock thickness for the viscous problem, 
T = Deoshe Vii (5.18) 


where T is the thickness, defined as the spatial distance from 
trough to peak of the waveform, and A = 2(1 +0)/mr. As an 
example, consider the case when [ = 100 and o = 1; then T = 
aM cosh +51 © -UG. ime this case, ne changes from a maximum 
negative value to a maximum poSitive value over a change in x 
Si. 0c. Since eat Y eee the change in u over this distance 
1S approximately twice the Mach number, which as noted before 
is taken to be of oon Then an approximate value of ae 


may be taken as Au/Ax x.25; therefore, while the I, term does 


become important, it does not dominate. For larger values of 
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ieee wili become increasingly important, will in fact domi- 
nate the equation, and will require smaller step size near a 
Beneck to correctly estimate eos 

Another consideration in the finite difference scheme is 
the handling of boundary values. Treated as an initial-value 
problem, there would be no problem in starting the solution, 
Since € and Ee would be everywhere prescribed for t = O. In 
mimeo case, Only one boundary condition is available for x = O, 
Piemoecher condition being that € > O as x + », which is of no 
import here except that it reguired the ruling out of solutions 
which grow in space. Because two "columns" of information are 
needed at the boundary to start the solution for g, another 
condition is required. Recalling the discussion in the previ- 
ous chapter concerning the modelling of boundary conditions, 
aimempcoblem of requiring an extra condition could lead one to 
expend great effort in attempting to develop a consistent and 
meaningful extra condition, when in fact if one simply applies 
the method of parametric differentiation as it is defined, the 
second condition will become apparent. The point is that the 
boundary is a real thing, not some abstraction in parameter 
Space, and consequently the source excitation has no depen- 
dence on w whatsoever; therefore one can in fact, without wor- 
rying about what Ey 1s at the boundary (which, incidentally, 
would overspecify the problem), write another boundary condi- 


iron , 
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g, (0, i) =O (5.19) 


This completes the information needed to march out the solu- 
tion of the problem. It is interesting to note that because g 
and g, are zero at the boundary, the only non-zero term in the 
finite difference cube on its first application is the inhomo- 
geneous term. Without this term, there would be no solution 
because all terms in the finite difference cube would be zero. 
Again, this iS consistent with intuitive thoughts about the 
form OL the solution for g; it should in fact be very small 
near the boundary so that when the quadrature is performed to 
calculate the particle displacement, the net change in the 


value of & from the inviscid to the viscous solution is small. 
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VI. RESULTS 

The practical realization of the preceding development is 
a computer program which solves the non-linear plane wave pro- 
pagation problem in a viscous fluid and provides output which 
allows one to readily compare the results obtained with this 
method, and those obtained by Keck and Beyer [4]. Their solu- 
tion was chosen as the point for comparison because it was 
Simple and easy to calculate, and because this analysis had to 
be conducted on a shoestring budget. Unfortunately, this did 
not allow good comparisons for the higher ranges of I because 
the Keck/Beyer solution is not sufficiently accurate [6], 
dacking harmonic content beyond the sixth harmonic. One can 
argue though that the results for large T appear to be consis- 
tent with one's expectations, and that the method appears to 
generate a plausible progression of results as the value of I 
is decreased (corresponding to an increase in WwW). An annotated 
version of this computer program appears in Appendix A. It is 
worthwhile to note at this point that because of the hyperbolic 
nature of the problem, the domain of dependence for a given 
point becomes larger and larger the further away from the 
boundary that the point is located. The result is that if one 
desires a solution value N discrete steps away from the boun- 
feany, N- + 2N + 1 points are required in the solution grid. 
This is a fundamental limitation of this method as it is pres- 


ently formulated, and may be difficult to improve upon. 
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In this analysis the Mach number was taken to be rather 
large, €« = O(.01), in order to minimize N while still being 
able to generate useful results. B/A was taken as 6, which is 
the approximate value of that quantity for water; thus £6 = 4. 
MemmemaGiven 8 product, the discontinuity distance, X, will 
occur at a point which is 1/278e wavelengths from the boundary, 
regardless of frequency. Thus, for Be = .04, X ~ 4A. In water, 
a Mach number of .01 corresponds to an acoustic pressure given 
by the expression p = PoCoe or p ~*~ 225 atmospheres, or in deci- 
bel units, SPL = 267 dB re 1 uPa, which is guite a strong sig- 
nal. The values of wy have been varied from zero to .01 which 
would reflect a change from inviscid to highly viscous propa- 
gation. Recalling the form y = Vb/Cod, it is clear that by 
ming YW, the frequency for which the solution is valid is al- 
so fixed for a given fluid. To put these results in a proper 
perspective, a value of J = .01 in water would correspond to a 
frequency of approximately five gigahertz. However, the re- 
sults for smaller values of wt represent more reasonable fre- 
quencies. Another interpretation of the results for y = .01 
is that a solution has been obtained for the fluid with an ex- 
tremely high viscosity coefficient such as glycerin, a solu- 
tion which has not previously been obtained analytically, be- 
cause of the restriction a/K <<1, which is necessary to re- 
duce the equations to analytically tractable forms. One might 


conclude from these results that there is nothing significantly 
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@ifferent in the form of the solution for such a liguid. In 
any case, the results which follow will indicate the potential 
of this method for application over a great range of values of 
the non-linear and viscous parameters. 

Prgunres 5a, b, C, and d are plots of the particle dis- 
placement calculated by parametric differentiation, the par- 
ticle velocity calculated from this solution, the difference 
between the inviscid (Fubini) and the viscous (Keck-Beyer) 
displacement solutions, and the difference between the invis- 
Smaeand viscous (parametric differentiation - P.D.) displace- 
ment solutions, respectively, for e = .01, and » = .01 (T = 
1.27). The two plots of the calculated difference are the ba- 
Sic results used for evaluation of the parametric differentia- 
tion solution. The differences are calculated according to 
Equation 4.16, and one would expect the calculated differences 
to be in phase with the displacement solutions; since viscos- 
ity should reduce the absolute amplitudes, the difference 
should be maximum when the Fubini solution has a maximum and 
minimum when the Fubini solution has a minimum. The zeroes 
will not necessarily coincide since the effect of the non-lin- 
earity is to lengthen the compressive portion of the waveform 
while reducing its magnitude and to shorten the rarefactive 
portion while increasing its absolute magnitude; viscosity will 
tend to keep the waveform sinusoidal, thus maintaining the 


length of the two portions of the waveform. Thus one could 
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hardly expect the zeroes to coincide. 

Since [T = 1.27, one can infer that the Keck-Beyer solu- 
tion should be accurate, and also that some harmonic content 
emould be evident in the particle velocity solution, assuming 
rT = 11s an accurate representation of the incipience of im- 
Bertant non-linear effects. In analyzing Figures 5c and 5d, 
One iS immediately struck by the fact that although the calcu- 
lated differences seem to be of the same form, there seems to 
be some sort of ramp type behavior entering into the P.D. dis- 
placement solution. If one were to try to place a physical 
interpretation on this result, it would be that the ability of 
the medium to restore itself elastically has been exceeded due 
to the large displacement amplitudes and that therefore a per- 
manent displacement has been induced, just as if a spring had 
been stretched beyond its elastic limit inducing a permanent 
strain. Due to the fluid nature of the medium and because no 
analytic solution predicts such behavior, such a physical in- 
terpretation may not have much meaning. An attempt will be 
made to offer a possible explanation for this result ina la- 
ter paragraph, but this remains an unsolved problem in this 
work. The average numerical loss in the positive portion of 
the waveform is clearly greater than that of the negative por- 
tion, rather than being equal as one would expect. 

The velocity solution, on the other hand, seems to be 


quite normal. One can observe that though non-linear effects 
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are producing some steepening of the profile, viscous dissipa- 
tion does in fact reduce the amplitudes properly; the ampli- 
tude of the last peak is 3.1 dB down vice 3.3 dB down by the 
Keck-Beyer solution, a difference of about 6%. The last point 
in the velocity profile 1s numerically inaccurate because of 
the lack of sufficient points at the tip of the solution mesh. 
The small “bump” in the velocity profile is not easily ex- 
mained, but probably arises from numerical inaccuracies. It 
seems remarkable that with little direct effort aimed at re- 
ducing computational inaccuracies in the numerical scheme, the 
difference in the solutions is but 6%, remembering all the 
while that this difference may not represent an error at all, 
but rather a more accurate solution produced by a method which 
does not resort to common approximations. One may wonder at 
this point why an equation in particle displacement vice one 
in particle velocity was treated, the answer being that the 
equation in particle velocity would still have required one to 
know the displacement solution in order to arrive at the vari- 
able coefficients which would have appeared in the parameter 
space equation, and thus would probably have been no essential 
Simplification. Another aspect of the P.D. solution which is 
worthy of note is the steepness of the velocity profile. The 
Keck-Beyer solution for this value of [ looks essentially like 
a decaying sinusoid, while the P.D. solution appears to indi- 


Cate that the non-linearity remains fairly important, result- 
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ing in some steepening of the profile. This suggests that 
meehaps | = 1 may not be a sufficiently strong indicator of 
Eeewancipience of non-linear effects. 

Results analogous to those of Figure 5 have been obtained 
for different values of the problem parameters, but with the 
number of points utilized in the calculation maintained at 
MeeeeolL. Figures 6 and 7 are the plotted results for ¢€ = .0l 
como = .0O01 (F = 12.7), and for ¢ = .01 and wy = .0000463 (T 
= 275) respectively. Figures 8, 9, and 10 are the results for 
meee O2 and ) = .01 (TF = 2.55), w = .001 (TF = 25.5), and wb = 
-0000463 (fF = 550) respectively. In addition, calculations 
have been conducted for ¢€ = .0075 and ¢€« = .005, but are not 
shown here. These results are not much different in form from 
those discussed previously, but they uncover some problems with 
the solution technique as well as help demonstrate that the 
method is in fact a potentially useful one for investigation 
of non-linear wave propagation problems, since it successfully 
solves the equation of motion in a predictable and consistent 
manner for various values of [. The range of the Be product 
had to be restricted as it was because of the need to minimize 
Sempucer costs and thus the number of points at which the dif- 
ference equation was to be solved. 

mr ting attention to Figure 6, first recall that with [I 
fmeeecne Keck-Beyer solution is probably only marginally ac- 


Curate. The displacement solution again exhibits the ramp type 
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behavior, although it is only readily clear when analyzing 
Figure 6d, which shows a greater average loss in the positive 
portion of the waveform than in the negative portion. Figure 
6c Shows that the Keck-Beyer solution is in fact no longer ac- 
curate near the discontinuity distance. The velocity waveform 
appears to be consistent with expectations with only slight 
attenuation because of the relative smallness of the attenua- 
mom coefficient over the spatial length of the calculations. 
However, small signal attenuation, l1.e., that attenuation which 
would occur if no harmonics were non-linearly produced, would 
mecaact a local Mach number of .0092 at the discontinuity while 
in this case a value of .0096 was obtained. This result is 
not consistent; one would expect greater attenuation in the 
non-linear case because the generation of harmonics will re- 
duce the magnitude of the primary in addition to small signal 
attentuation, while the harmonics themselves will be attenu- 
ated even more rapidly because of greater attenuation coeffi- 
cients. Unfortunately, in this case there is no easy explana- 
tion of the result; however, it did prompt some thoughts about 
what the Fubini solution looks like and what the changing 
Shape of the waveform should imply with respect to the instan- 
taneous values of particle velocity near the discontinuity 
@aestance. 

From the point of view of energy conservation, one would 


expect that in a fluid without viscous dissipation, the time 
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Or space average energy density in a plane wave train propagating 
isentropically will remain constant. The implication of this 
is that as the finite amplitude waveform steepens and eventu- 
ally assumes a sawtooth shape, the peak values of the particle 
velocity must increase in order to maintain constant average 
energy density in the wave train. For a sawtooth wave to have 
the same energy density as a sinusoid of unit amplitude and 
identical period, its amplitude must be approximately 1.22. 
Realizing that a fully formed sawtooth develops at a distance 
/Mmeom the boundary of about 3.5xX [3], it would still be expected 
that the peak value of the particle velocity at the discontin- 
uity distance would be greater than that at the boundary. The 
MerseessOlution (Fubini) calculated in this study is accurate to 
five places in displacement and three in velocity, and the 
peak particle velocity at the discontinuity distance is ident- 
meee tO that at the boundary, implying that there in fact has 
been an energy loss in an inviscid solution, and it is certain- 
ly not clear where this energy has gone. In any case, for a 
Meete'z./, non-linear effects are definitely important, and it 
is entirely conceivable that the P.D. solution has somehow 
recognized an error in the base solution and has yielded a 
quadrature value which correctly accounts for the steepening 
Semene profile and thus predicts a particle velocity which is 
greater than that which would be obtained with a sinusoid and 


Stemeesignal attenuation. This explanation is of course sup- 
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position since there are too few points calculated to obtain 
a good space average energy density to verify an energy loss 
which is greater than that predicted by small signal attenua- 
tion. It indicates an unsolved problem and at the same time 
Suggests that parametric differentiation predicts the proper 
solution to the problem. Blackstock [9] suggests the mechan- 
ism for dissipation in an inviscid fluid after a shock has 
formed; here, there 1S no reason for such dissipation since 
Pietenare no shocks in the first zone of non-linear propaga- 
iON . 

Figure 7 shows much the same results as those of Figures 
meena 6. In this case, [ = 275, and the Keck-Beyer solution 
is clearly not converging to the proper result as evidenced by 
Figure 7c. Again the ramp behavior is clear in Figure 7d, 
though the magnitudes of the difference values are understand- 
ably smaller. The peak value of the calculated Mach number at 
the discontinuity distance is now .0196 which is in fact great- 
er than the peak at the boundary. The differences between the 
inviscid and viscous displacement solutions are now so small 
that on the scale of these plots, the viscous solution appears 
mdemtacal to the inviscid solution. This is not remarkable, 
but it is noteworthy that the plotted differences, though 
small, produce a smooth curve rather than one that jumps from 
point to point. Figures 8, 9, and 10, which represent the 


calculations for a source Mach number of .02, show results 
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which are similar to those shown in Figures 5, 6, and 7, and 
do not need to be analyzed with as much detail. However, in 
Pagure 8, with [T = 2.55, the peak Mach value near the discon- 
Mmaenuaity distance is .0147 vice .0137 according to the Keck- 
Beyer solution for a difference of 7%, which is consistent 
with the 6% difference obtained in the case ¢€ = .01, fT = 1.27. 
It is also consistent with the idea that peak Mach numbers 
should increase in the non-linear case since small signal at- 
tenuation would predict a Mach number of .014 at the same 
point for a purely sinusoidal wave. The ramp behavior is evi- 
Meme in each of the Figures 8d, 9d, and 10d.Finally, these results 
also provide a counter example for the idea that the peak 
Mach number should increase near the discontinuity distance in 
mmeecase Of strong non-linearity, 1.e. when [ is greater than 
Mee the particle velocity plots of Figures 9b and 10b indi- 
cate peak Mach values which remain less than the peak at the 
boundary; in Figure 9b the value is .0193 and in Figure 10b it 
is .01998. Small signal attenuation would predict a value of 
Meee £Or the result of Figure 9b. Energy considerations dic- 
tate an increase in the peak Mach values and in this case it 
does not occur. However, the spatial step size used to arrive 
at the plots of Figures 6 and 7 was .04, while in the calcula- 
tions used to obtain the plots of Figures 9 and 10, the spa- 
tial step size was .02. Recalling the finite difference ap- 


proximations utilized, this difference in step size may be the 
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explanation for the difference in results. The truncation er- 
ror for the second-order difference approximations is of ox?) 
Hieretore, it is certainly conceivable that the step size of 
504 is not accurate enough to yield good results for particle 
velocity. At the same time the considerations of energy can- 
not be ignored. It is likely that the increase in the peak 
Mach number in the results of Figures 6 and 7 is a result of 
mamerical inaccuracy, but this is certainly not entirely clear. 
The other calculations conducted with Mach numbers of .0075 
and .005 produced results which were consistent with the other 
results obtained, but they are not presented here because the 
spatial step size was sufficiently large that the results are 
somewhat questionable because of truncation error. 

There remains a question with respect to the curious ramp 
behavior in the particle displacement solutions. In examining 
the plots of the difference between the inviscid and viscous 
solutions, this behavior is always clear, but would not neces- 
Sarily be clear at all in analyzing the particle displacement 
results, particularly as ) becomes smaller and smaller. vw is 
a direct measure of the viscosity of the assumed fluid rather 
than a relative measure of viscosity such as [. Recalling that 
aw of .01 would imply a highly viscous liquid or a very high 
frequency wave and that the approximation upon which previous 
authors' results are based is a/K << 1, which implies weak 


viscosity, the reason for the ramp may become clear. Previous 
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analyses, in making this approximation, may have neglected a 
rather interesting physical interpretation of the non-linear 
loss mechanism. For small wt, the values of the calculated 
differences are so small compared to the values of particle 
displacement that this effect would not be noted when solv- 
ing for displacement or velocity without consideration of the 
differences between the viscous and inviscid solutions. In 


fact, in making the assumption a/K << 1, the values of jy are 


apparently restricted to relatively small values; a/K = wee 
O 
= me ppand Y= a8. , therefore it is clear that previous 
O O 


analyses have not considered values of on the order of those 
considered here. For large W and thus highly viscous propaga- 
tion, previous results are not valid, and therefore if the 
ramp behavior has a physical explanation, it could not have 
been noted before. 

In the parameter space Equations 4.6 or 4.9, there are 
Peueecerms which involve the dependent variable, Ce 


g and Sy.) which result from similar terms in the displace- 


One 
ment space Equation 3.9. Gee 1S a term which relates to the 


elastic properties of the medium, to the inertial proper- 


Jet 
ties, and onan to the viscous properties. The Sy. term implies 
that there exists some modification of the elastic properties 
of the medium. Recalling the discussion in the previous sec- 


tion concerning the effect of this term in the portions of the 


waveform where the velocity gradient u, assumes large values, 
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the source of the ramp behavior may become evident. When oo 
isetarge, the Sy. term will dominate the equation. In the in- 
Peecerd SOlution, oy beeoMNcoulMrimtte at the first Shock. 
Semee the inviscid solution is in fact the starting solution 
and the quadrature from parameter space to displacement space 
is accomplished with an integration with respect to yp, it is 
obviously not unreasonable to expect oe to assume large val- 
ues when is small; also, one would expect a to assume nu- 
merical values which are more accurate for the current y step 
of the calculation if the spatial step size was small in the 
regions of large Ser In any case, if the Jy. term dominates 
maeweauation in a specific region, the g solution in that re- 
gion will vary linearly with x, hence the € solution will vary 
in a similar manner, and hence the ramp. 

mrethe absence of viscosity, it is clear that the fluid 
Mearercles oscillate about an equilibrium position; with vis- 
cosity, realizing that the g,, term 1s one which, in essence, 
alters the elastic constant of the system, it 1s postulated 
that the following phenomenom occurs. Because of the large 
velocity amplitudes, the phase velocity increases in the com- 
pressive portion of the wave and decreases in the rarefactive 
portion; this 1s well known. These large particle velocities 
produce correspondingly large particle displacements, but now 
Mem rluid particles, instead of being purely elastic, do not 


have all of their kinetic energy at zero displacement converted 
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to potential energy at full displacement. Some is lost to 
Small signal attenuation of both fundamental and harmonics, 
and the rest produces an irreversible displacement of the 
fluid particles; that is, work is expended on the fluid system 
to produce a displacement of the equilibrium positions of the 
mana particles. Unfortunately, this effect can only be pos- 
tulated here because a precise energy balance has not been 
Performed. In order to calculate the specific energy lost in 
producing a given displacement at a given point, the net ac- 
celeration on the fluid particles which produced this displace- 
ment would be required, and this calculation has not yet been 
performed. 

The reason for this apparent displacement of the equili- 
brium positions of the fluid particles may be explained in a 
more rhetorical manner as an interaction between non-linear 
and viscous effects. The non-linear effect is to alter the 
elastic constant thus changing the speed of propagation of the 
wave, resulting in the steepening of the waveform. Now, be- 
cause the oncoming compressive portion of a given cycle of the 
wave is travelling faster than the preceding rarefactive por- 
tion and because viscosity is acting to slow the velocity of 
Permit iuid particles in the rarefactive portion even more, the 
Meamercles, after passage of the given cycle, do not return to 
their former equilibrium position, but are instead minutely 


displaced. The oncoming compressive portion of the wave is 
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'filling’' the preceding rarefactive portion and consequently 
the fluid particles in the rarefactive portion do not have the 
Opportunity to recover their initial equilibrium positions. 
Such an effect would obviously tend to create a vacuum behind 
the wave train, but this can be justified by requiring that 
fluid flow in from outside the region excited by the source. 
Finally, there are some matters concerning the numerical 
method which require attention. The stability criterion de- 
rived in the preceding chapter was arrived at in a rather un- 
conventional manner. Some numerical experimentation was per- 


formed to determine just how sensitive to the ratio of step 


sizes the solution was. The exact stability criterion was 
fae Nt / (1 + 2e)-P Go 
_ Ox 
Pemeneimplies for d&/dx positive, e¢ = .01, and B = 4 that Ax 


fmeey7). 08 near the first shock. This should be an absolute 
limit on the spatial step size since 9&/9x will have its maxi- 
mum value at this point. The equations were solved for Ax=.04 
and At=.04, and the numerical solution became highly unstable 
beyond the second cycle of the waveform. Then At was set equal 
to .045 which would be just slightly greater than the At re- 
quired by the stability criterion for Ax=.04, and the solution 
was stable throughout. All subsequent calculations were con- 


ducted with a convenient At selected so that it was large enough 





aoe 


Bem the given Ax to satisfy the stability criterion. 

Another matter of concern in the calculations is the ef- 
meeteror the step sizes, in x and t and in ww as well, on the 
meecuracy Of the solution. In this analysis, the base solution 
wes Calculated accurate to five significant figures, obviously 
mmeeing the accuracy of the calculated solution to the same 


extent. The base solution required 51 harmonics at the dis- 


Semeinuity distance for this accuracy. The spatial step sizes 
used for the various solutions were Ax = .08, .06, .04, and 
-02. Since the truncation errors in the second order differ- 


ences are of o( Ax) 7), 1t would be reasonable to assume that the 
Paweulated values of g would be accurate to two or three places 
depending on the step size used. However, to return to dis- 
placement space, the values of g are multiplied by a difference 
between two values of w, and hence for a maximum w of .01, if 

g is accurate to three significant figures, the differences 
calculated in the quadrature will leave the displacement solu- 
lmom accurate to five significant figures, neglecting errors 
introduced by the discretization of the solutions with respect 
to WwW. This then would be consistent with a base solution cal- 
eemeaeea tO four significant figures, for Ax=.02 and marginally 
consistent for Ax=.04. This is the reason why the results for 
€ = .005 and « = .0075 were not presented; their accuracy could 
Only be considered marginal at best for the maximum yj value. 


These were the solutions which employed step sizes of .06 and 
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moo. Roundoff error was not considered to be of significant 
concern in this analysis with machine accuracy to seven Signi- 
ficant figures. The solutions also did not seem to be sensi- 
tive to the step size in time as long as the time step was 
immege enough to satisfy the stability criterion at all points. 
Finally, the step size in ~ seemed to matter very little in 
the form of the results if not in the exact numerical results. 
memes Obvious that in integrating across orders of magnitude, 
the integration must be handled carefully, but experience dem- 
Onstrated that values of g at a given point changed little as 
W was varied, and thus the error introduced by utilizing a 
relatively large step size in wW should be small. Consequently, 
the following scheme was adopted: 
1) the first ~ value was the limiting value Ye 
2) the second wW value was that which when used in 
the quadrature yielded difference values which 
were just within machine roundoff accuracy 
3) each subsequent step in the order of magnitude 
of ~w was broken into an algebraically increasing 


number of steps on a logarithmic basis. 


For example, with Bo = O, Vy = momo there would be two steps 


= 5 5 


between ) = 10° and w = 10 ~, with (ee oe 10 ~ and WV. 


ze non” , Enuec steps between i = Ona eercl yp = non) aes 


oq mon, Ve = .463 x cme and Ve = 104, etc. There is no 
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clear method to determine the proper step size in w, but this 
one produced consistent results. In general, it would be in- 
meerOopriate 1f g varied in a highly non-linear manner with }. 
The final comments concerning the numerical results are 
made with respect to the formulation of the problem. Equations 
4.6 and 4.9 were the parameter space equations developed from 
the equation treated by Keck and Beyer and from the parent e- 
M@eaelon without approximation, respectively. Early calcula- 
tions were performed with Equation 4.6 and later ones with E- 
quation 4.9. For ¢ = .01 and wy = .01, the maximum percentage 
difference between the calculated solutions by Equations 4.6 
and 4.9 was approximately 4%, demonstrating that the equation 
solved in previous work does in fact produce very little error 
for Mach numbers up to .01. Finally, the two forms of the 
boundary values for g appeared to make little difference in 
the final result; producing a maximum percentage difference in 
the calculated values of particle displacement of approximate- 
[ieeico between the case utilizing the ‘natural' boundary condi- 
tions for g and the case in which artificial dependence on jp 


was introduced at the boundary. 
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VII. CONCLUSIONS AND RECOMMENDATIONS 

A theoretical study of non-linear plane acoustic wave 
propagation in a viscous fluid has been conducted utilizing a 
little known mathematical technique, parametric differentiation, 
in Order to demonstrate its utility when applied to acoustic 
propagation problems, and to obtain whatever physical insight 
into the problem the results may present. It has been demon- 
meraced that the method can in fact yield results which are 
consistent with other analytic solutions to the problem for 
certain arbitrarily selected values of the viscous and non- 
linear parameters. No restriction was required in the selec- 
tion of these values by the mathematical nature of the method, 
and there is no reason to believe that the method could not be 
applied for any values of the parameters. Restrictions which 
may exist in the application of the method are computational 
ones, specifically: 1) the results can only be as accurate as 
the solution to the problem for yp = Voi 2) ew aceiuiacy may be 
further limited by the finite difference approximations util- 
mea; 3) accuracy may be limited by the ability of the computer 
to solve the finite difference scheme for a large number of 
points with a given time and memory budget; for instance, dou- 
ble precision will require twice the computer memory; and 4) 
the hyperbolic nature of the problem requires calculations to 
be conducted throughout the domain of dependence of a given 


Meme tO arrive at a solution at that point. These restric- 
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tions may render the method unusable for calculations involv- 
igelanlge propagation distances, but they do not detract from 
its ability to provide physical insight. Further research may 
provide a reasonable technique for application of the method 
when large distances are involved. The solution of the prob- 
lem by parametric differentiation has allowed the elimination 
of the restriction a/k << 1 heretofore assumed in typical pro- 
meeation Studies. Elimination of this restriction has produced 
results which are valid for highly viscous fluids as well as 
demonstrated that the results for such fluids are essentially 
Similar to those for less viscous fluids. It has allowed the 
formulation of a postulated behavior for all viscous fluids, 
which may explain physically why there are non-linearly in- 
duced energy losses above and beyond those associated with 
viscous attenuation. Specifically the solutions produced by 
this method predict that some of the energy contained in a 
finite amplitude wave in a viscous medium iS expended in pro- 
ducing a displacement of the equilibrium positions of the fluid 
particles. These results have also been obtained cheaply; for 
101 spatial points, the computer time required to produce a 
result for single values of wW and the fe product is approxi- 
mately .15 minutes with 400 kilobytes of computer memory. 

This is not to suggest that this work proves beyond a 
doubt that the method is a useful one. Consideration of prop- 


agation once a shock has actually formed has not been attempted. 
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If the shocks were in fact discontinuities as predicted by the 
inviscid theory, the technique probably would not be useful 
beyond the first zone of propagation. However, since the wave- 
ferm 1s in fact continuous through a shock of finite thickness, 
it is reasonable to expect that with denser sampling of the 
waveform through the shocks, good numerical results could be 
Seecained. The next step in a continuing study of the solution 
of the problem through parametric differentiation would be to 
Select Blackstock's weak shock solution [3], which connects 
the Fubini solution [1] with the Fay solution [2] in its in- 
esciad limit to yield a solution valid to the end of the sec- 
ond zone of propagation, as the base solution for ve = O, and 
to attempt to extend the solution into the region of periodic 
shocks. With this result a precise energy balance calculation 
would demonstrate whether or not energy is lost in producing a 
permanent displacement of the fluid particles. 

More fundamentally, it would be useful if a mathematical 
analysis of the technique were to be conducted in order to de- 
termine why it apparently correctly predicts solutions for many 
different classes of problems, and what if any are the theor- 
etical limitations to its application. Beyond this, there are 
more interesting problems in non-linear wave propagation to 
which the method could be applied. The first, and one which 
can be readily adapted to the computer program already written, 


would be an investigation of the solution of the Burgers equa- 
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tion as IT is varied. Having noted previously that this solu- 

tion is difficult to calculate as TF becomes much greater than 

memor for | > 50 because of slow convergence, it is suggested 
efefertne analytic solution be calculated for f = 1. Then with 
the Burgers equation in the form solved by Blackstock [5], 


with V = ee 


vee- VV = / V_, Vio, = sinw eis 
a y ee) = Savy ol) 


differentiate with respect to [ to obtain, 


ceo 9 = 1 Vy glo, y) = 0 (7e2) 
where o 1S a spacelike, and y a timelike variable. This equa- 
tion would be perfectly suited to solution on the mesh used in 
this analysis because there would be no value at the itl1,j-l 
point, and the scheme would be completely explicit. The solu- 
tion for Tf = 1 could then be extended to arbitrary values of [ with- 
out significant computational problems. Because of the space 
scale of Equation 7.1, this technique will allow easy calcula- 
tion of the solution in both the first and second zones of 
propagation. 

Finally, there are the two problems of finite amplitude 
wave propagation which are yet to be solved analytically, those 


Of propagation of cylindrical and spherical waves in viscous 
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fluids; parametric differentiation may provide an invaluable 
tool in solving these problems. The fundamental equations of 
motion do not lend themselves to analytic solution, there be- 
ing no simple transform such as that used by Blackstock in the 
plane wave problem to reduce the equation of motion to a Bur- 
gers equation [18]. Indeed, the inviscid forms of the equa- 
tion have no exact analytic solutions which have been published. 
Blackstock [9] has obtained approximate solutions to the cylin- 
O@rical and spherical wave problems for an inviscid fluid, one- 
dimensional wave propagation, and r >> 1, which are very sim- 
ilar in form to the Fubini solution; these results were ob- 
tained through manipulations of the equation of motion which 
yielded an inviscid Burgers equation. Such solutions could be 
used as one base solution for a solution to the viscous prob- 
lem by parametric differentiation, realizing that they are on- 
ly approximate. Other authors [12][19]} have produced particu- 
Ma~esOolutions for spherical and cylindrical wave propagation, 
but with an accurate base solution, parametric differentiation 
could be used to generate many solutions, for a fixed Be pro- 
Guct, with relative economy. Consequently, the apparent course 
Semoeudy tO pursue for non-linear cylindrical and spherical 
wave propagation is the development of improved base solutions 
Or Vo = O and the application of parametric differentiation 

to these solutions, producing a large number of solutions which 


may then allow for some generalization of the results. Such a 
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research program could lead to a significant gain in the under- 
standing of the physical mechanisms which are collectively des- 


mee as non-linear attenuation in viscous fluids. 
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Figure 8c. Difference Between Fubini and Keck-Beyer 
Displacement Profiles, Be = .08, I’ = 2.55 
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